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We present the results of a Molecular Dynamics computer simulations of a two component isotope 

• mixture of Lennard- Jones particles, monodisperse in size but different in masses, at a fixed average 

' density and temperature. We study changes in properties that result from mass heterogeneity, by 

, measuring the pair distribution function, difi^usion coefficient, velocity autocorrelation function, non- 

. Gaussian and isotope effect parameters, as functions of the degree of mass difference. Our results 

, ' show that if static properties are not influenced by a variation of mass variation, the dynamic 

I properties are significantly affected and even exhibit the presence of "critical" values in the mass 

, . difference. We also demonstrate that our model gives a simple contra example to a recently proposed 

' a universal scaling law for atomic diffusion in condensed matter [M. Dzugudov, Nature, 381, 137 

: (1996)]. 
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INTRODUCTION 



The aim of this paper is to study, using Molecular Dynamics (MD) computer simulations, the 
effect of having a mass difference in particles on properties of liquids. The influence of mass 
polydispersity on dynamical properties when simulating the dynamics of realistic system having 
. size polydispersity was recently demonstrated by Poole and one of the authors (N.K.) Computer 

simulation methods, MD in particular, gives unique possibility to isolate and study the role of a 
single parameter for systems of any complexity that is much more difficult to accomplish (or even 
'"^ ' impossible) in an actual experiment. 

^ I Here we consider a model that allows us to isolate completely the effect of only changing the 

O ■ masses of the particles by keeping all other characteristics of system unchanged. Our model is the 

simplest possible isotope system: a two-component mixture of isotopes, particles of the same size 
but different masses, consisting of 50% particles A with mass uia = mo — Am and 50% particles B 
with mass ms = mo + Am, so that no matter how we change Am, the average mass of particles, 
^ j , total mass and average mass and number densities of a system remain the same. In this sense, our 

■ model differs from these studied previously. For example, Ebbsjo et al. in and Bearman and 

Jolly in studied isotope systems in which one species had a fixed mass and the mass of another 
species varied. Consequently, this affected the average mass and mass density. Other authors 
considered the diffusion of solute particles in the limit of infinite dilution in a solvent. 
We perform the equilibrium MD simulation of a 3-dimensional system of iV = 4000 particles 
interacting via the shifted-force Lennard- Jones (LJ) potential a modification of the standard 
LJ potential 



V{r) = As 



cr\12 /crN6 
r / \r 



(1) 



where e characterizes the strength of the pair interaction and a describes the particle diameter. 
Both e and a are constant for all particle pairs. In the shifted-force LJ interaction, the LJ potential 
and force are modified so as to go to zero continuously at r = 2.5(7, and interactions beyond 2.5ct 
are ignored. We chose the Lennard-Jones potential because it is very popular in MD studies and 
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has proved to be a good approximation, for example, for liquid argon. In addition, there is a large 
amount of data, both theoretical and experimental, available for comparison [lUi fill Ha. 

In our work we use reduced units. Energy is expressed in units of e, mass m in units of mo, 
length in units of tr, the number density of particles p in units of cr^"^, and temperature T in units 
of e/k, where k is Boltzmann's constant. Time t is expressed in units of a/ mocr^/e. In these units 
the time step used for integrating the particle equations of motion is 0.01. The initial distribution 
of velocities was Maxwellian. 

After equilibration, all quantities are evaluated in the microcanonical ensemble (a constant 
NVE). We present data for p — 0.75 and T = 0.66 that is not far from a triple point of one 
component fluid (pt = 0.85, Tt = 0.76 14]) where the cooperative motion of liquid particles is more 
prominent. We conduct simulations for different Am from to 0.7 with the step 0.1. 



II. PAIR DISTRIBUTION FUNCTION, DIFFUSION COEFFICIENT, AND 
VELOCITY AUTOCORRELATION FUNCTION 

The pair distribution function g{r) that characterizes the average liquid structure [l5l | is shown 
in Fig.^ We calculate g{r) for systems with different values of Am and for each species A and B of 
these systems. We find that the pair distribution function g{r) is the same in all cases, confirming 
that the static properties of the system do not depend on Am jl^. 

Fig-Hlpresents the dependence of {\r{t) — r(O)p) (the mean square displacement of a particle) on 
t for different Am. As Am becomes larger, (|r(<) — r(0)|^) has a steeper slope. In addition, in Fig. 
I^lwe plot three graphs that show the mean square displacements of particles A, B, and the total 
(|r(i) — r(0)|'^) for Am = 0.7, as an example. The curves for particles with mass itla — ttlq — Am 
lie above and for particles with mass ms = mo + Am lie below the curve for the total mean square 
displacement. This demonstrates the fact that lighter particles not only move faster between 
collisions, but also diffuse faster than heavier particles. It is interesting that for small t the graph 
of the mean square displacement of a lighter particle has a steep slope which decreases and becomes 
constant for t > 0.15, whereas the graph of the mean square displacement of a heavier particle 
increases for < t < 0.15. 

We calculate the diffusion coefficient D in two ways: using the Einstein relation 

. 'I--'" -;'»>!'' (2) 

and using the Kubo formula [l5| 



DK = \j^ (v(0) • v(0) di, (3) 

where (v(0) • v(t)) is the velocity autocorrelation function. Here and below subscripts E and K 
mean that D calculated using the Einstein formula and using the Kubo formula. 

The results are presented in the Table Da(b) is the diffusion coefficient in scaled units for 
species with mass m^(5), and D is the average diffusion coefficient of the mixture. At fixed p and 
T the diffusion coefficients D, Da, and Db increase as Am increases. 

De and Dk are equal within the computational uncertainties; however we find that calculations 
using Eq. (O are more accurate. The same conclusion was made in Ref. This is because Eq. 
© involves integration where negative and positive parts of (v(0) • v(t)) partially cancel out and 
because of uncertainties in the long-time behavior of the velocity autocorrelation function, whereas 
the graphs of (|r(t) — r(O)p) versus t are almost perfect straight lines for t > 0.15. 

From the Tabledwe can see that the relationship Da > D > Db always holds and the difference 
among Da, D, and Db increases with an increase of Am. This is shown in Fig. ^ The dependence 
of the diffusion coefficient versus Am is not linear. 
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Even though the diffusion coefficient does not change too much with a change of mass (if Am 
changes by 1%, D changes maximum by 5%), the velocity autocorrelation function varies signifi- 
cantly, both in shape and magnitude as a function of t. Fig. \^ shows the dependence on Am of 
the normalized velocity autocorrelation function 

(v(0) ■ v(t)) 

^^^^ = ^Fwr- 

It is interesting to observe the behavior of the negative part of The velocity autocorrelation 
function, basically, is the projection of a particle's velocity at time t on the initial velocity of that 
particle, averaged over all particles. If ip{t) becomes negative it means that a particle after a 
number of collisions, on average, reverses the direction of its motion. As we can see from the 
graph, all of the functions 4'{t) have a negative part but the minimum of il'{t) behaves differently. 
For Am = 0.1 and 0.2, the functions 4'{t) are very close to ipit) with Am = 0. As Am increases, the 
minimum shifts to earlier times. Its magnitude initially decreases, has a maximum at Am — 0.5, 
and starts to increase drastically for Am > 0.5. To reveal this difference we plot three curves on 
Fig. IHl which shows the velocity autocorrelation functions of particles A, B, and the total value 
for Am — 0.7. The curve for particles with mass m^ = mo — Am has a deeper minimum than the 
total and the curve for particles with mass m^ = mo + Am does not have a negative part at all. 
For a big difference in the masses m^ and niB , the heavier particles do not reverse their direction 
of motion on average, even in a dense fluid. 

III. NON-GAUSSIAN AND ISOTOPE EFFECT PARAMETERS 

The general non-Gaussian parameter Q;„(i) is defined as 0| 

where c„ = [1 • 3 • 5 • • • (2n -I- 1)]/3"' and (r^"(t)) are the ensemble average of the 2n*^ power of 
the particle displacements after a time t 



id ( r 

m 



1=1 



For systems in which motions of particles arc uncorrelated (for example, in an ideal gas) a„(t) = 
0. The deviation of a„(t) from zero serves to quantify the correlation of particle motions at 
intermediate time or, as it was shown in 0], the heterogeneity in masses for short times. 

Here, as in 0, we calculate a2{t): 

In Fig. |7|we plot a2{t) for different values of Am. For small Am the behavior of a2{t) is similar 
to that of a polydisperse system 0. For Am ^ 0, a2 does not start from 0. We observe the 
characteristic, intermediate-time peak of a2, at approximately t = 1, which increases in magnitude 
as Am increases. As in the combination of the non-zero, early-time behavior of Q!2, and 
the intermediate-time peak produces a minimum in a2 at approximately t — 0.12. However, for 
larger Am, starting from Am — 0.5, the minimum disappears, and a2(t) becomes a monotonically 
decreasing function. For these Am the values of a2{t) at i ^ even exceed the intermediate- 
time peak. It is interesting to observe how the total non-Gaussian parameter is related to ones of 



4 



particles A and B. Fig. |S1 shows three curves representing a2{t) for particles A, B, and the total 
when Am = 0.7. Both curves {A and B) start from (each component has particles of the same 
size). The curve for particles with mass = mo — Am has a maximum that is twice as high and 
is located at an earlier time than the curve for particles with mass niB ~ Tig + Am. The difference 
between these two curve increases when Am becomes larger. In addition, the curve for the lighter 
particle, after passing the maximum, completely coincides with the total a^it) and decays as 
according to [17| . The lighter particles are more mobile and give a larger contribution to the 
motion of the fluid. 

As we pointed out in 0, the non-Gaussian behavior of Q!2(/:) shows that the self-part of the van 
Hove correlation function Gs (r, i) 



Gs{r,t) = {^Y.^{r-\v,(t)~vm\)y (8) 

deviates from the Gaussian form. In the limit r, t — > 0, when particles move with a velocity 
Vi = Tj/t as if they were free, this corresponds to the condition that the distribution of velocities 
of particles is non-Maxwellian. The Maxwell-Boltzmann distribution for particle velocities is 




We can see that it depends on the mass of a particle. Even if for each species the distribution 
is Maxwellian the total distribution is not 



In Fig. Owe plot the distribution of velocities of all particles (solid line) and the Gaussian fit of 
this curve (dashed line). It is clear that this distribution deviates from the Gaussian form. It is 
narrower and has longer tails for large velocities. This is the contribution of particles with smaller 
mass. It is interesting that the same non-Gaussian behavior of Gs(r, t) was observed experimentally 
in ^3 for colloid suspensions. It could be accounted for the polydispersity in the mass of colloids 
used in that experiment (which was ~ 8%). 

Finally, we can use the formula for the non-Gaussian parameter q;2 at t ^ for the multicom- 
ponent system from 0: 



^2 = a2(i ^ 0) = — — n2-1' (11) 

where rii — Ni/N is the fraction of particles of species i. (Here, as in [J, the superscript "o" 
indicates the limit i ^ 0.) 

If we have a binary system with two species having masses mA — "nio ~ Am and tub — mo + Am, 
Eq. becomes 



'"■A I nB 2 

o _ (mo-Amy-' (mo+Am)'^ ^ _ inATlBX 

{l + {nA- nB)x) 



^?no— Am mp+Am^ 

where x ~ t^viilva^). In our case ua = = 1/2 and 



5 



al = x^. (13) 
We plot a2 versus Am in Fig. IIUI confirming this analysis. 

An additional parameter which could reveal the diffusion mechanism in isotope mixtures is the 
isotope effect parameter. It was shown by Parrinello et al. in fl^ that for dilute gases the ratio 
of the self-diffusion coefficients of a binary mixture is equal to the reciprocal of the square root of 
masses of the two species. However, computer simulations performed by Bearman and Jolly |^) 
show that the ratio of the self-diffusion coefficients varies with masses as 



^=(^\ (14) 

where 7 differs from 1/2 for higher densities and lower temperatures, when effects of collective 
motion is important. In Fig. 1111 we plot \ti{DaIDb) versus ln(r7iB/m^). The slope of the best 
fit line gives 7 — 0.081. This corresponds to well-known results (for example, therein Ref. pof'l. 
There is another parameter that serves as a quantitative measure of collectivity, the isotope effect, 
defined as [2fl 



E = ^d£M^. (15) 

Large isotope effects are interpreted as single-atom jumps via vacancies pl|. whereas small iso- 
tope effects indicate the presence of some collective processes. Kluge and Schober in [23| estimated 
the number of particles moving cooperatively as 



iV«i. (16) 

We find that for our binary mixtures the isotope effect E varies from 0.132 for Am = 0.1 to 
0.109 for Am = 0.7. It indicates an increase in the cooperative motion of particles with an increase 
in their mass differences. 



IV. CONCLUSION 



We have considered MD computer simulations of a toy model, a binary isotope mixture, that 
allows us to elucidate the effect of having a mass difference on properties of liquids. It has been 
shown that redistribution of mass (Am) among the two species, while keeping the rest of parameters 
fixed, does not affect static properties (the pair distribution function (7(r)) both for the system and 
separately for each species A and B and these are independent of Am). At the same time, all 
dynamic properties (the diffusion coefficient, velocity autocorrelation function, non-Gaussian and 
isotope effect parameters) are affected and exhibit Am dependence, which demonstrates the pure 
dynamical nature of mass heterogeneity. 

Moreover, the Am dependence exhibits some "critical" values such as in the behavior of the 
negative part of the normalized velocity autocorrelation function ■)/)(/:) (see Fig. which at Am = 
0.5 changes the initial tendency of decreasing magnitude of first minimum with increase of Am 
to the opposite, i.e. a drastic increase in V'(0- Such phenomenon, we believe, deserves further 
investigation and analysis. 

Many real systems have mass heterogeneities but, of course, not pure as in our model. Never- 
theless, even the hidden and probably not very strong effect of having a mass difference compared 
to other (size, interaction, etc.) has to create a pure dynamical contribution. As demonstrated 
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by our simple model, this pure mass effect does not affect the static properties of a system. In 
particular, this observation leads us to question the validity of some universal empirical formulae 
based on the conjecture that "atomic diffusion is an entirely geometrical phenomenon" The 
author of |22| claimed that the relationship between D and g{r) is governed by following equations 

D* ^ 0.0496-^% (17) 
where D* is the diffusion coefficient in the dimensionless form D* = DV^^a'"^, Te is given by 



- ia^g{a)p^nkT/m, (18) 

where cr is a position of first maximum of g(r) and ^2 is the two-particle approximation of the 
reduced excess entropy 

/>OC 

^2 = -27rp / {g{r) In [g{r)] - [g{r) - 1]} r^dr. (19) 
Jo 

According |2^, Eqs. H17I18I19|I are universal for equilibrium condensed atomic systems, both 
liquid and solid, regardless of the structures, interatomic interaction potential or the microscopic 
dynamical mechanisms involved and also valid for multicomponent systems with Te, o and 5*2 
corresponding to each type of constituent atoms. 

This universality contradicts results of our simulation because Eqs. (|17I18I19|) for the considered 
binary mixture give Eq. (|14|l with 7 — 0.5, whereas from Fig. ^2 the slope of the line of best fit 
gives 7 — 0.081. 

For our system the observed difference in the diffusion coefficients of two species is an entirely 
nongeometrical, nonstatic phenomenon. 
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Da(k) 
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D(K) 


1 
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5.4979 


5.4749 


5.4979 


5.4749 


5.4979 


5.4749 


0.9 


1.1 


5.5861 


5.4866 


5.4592 


5.7079 


5.5355 


5.6564 
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1.2 


5.()()77 


5.G;iG9 


5. 1710 
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5. .5 79 7 


5.18 12 


0.7 


1.3 


5.7393 


5.7635 


5.4802 


5.5991 


5.6108 


5.6169 


0.6 


1.4 


5.9216 


5.8336 


5.5586 


5.5582 


5.7045 


5.7818 


0.5 


1.5 


6.1117 


6.0612 


5.5852 


5.6850 


5.8550 


5.9024 


0.4 


1.6 


6.3923 


6.6003 


5.7286 


5.8319 


6.0542 


6.1483 


0.3 


1.7 


6.7438 


6.5358 


5.8589 


6.0439 


6.3819 


6.4645 



TABLE I: The dependence of the diffusion coefficient on masses of species. All D are multiplied by 10 
(in scaled units). 




FIG. 1: The pair distribution function g{r) of all species. 
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FIG. 2: The mean square displacement versus t for increasing Am (from the bottom to the top). 
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FIG. 3: The mean square displacement versus t for lighter, heavier particles and the total for Am = 0.7. 
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FIG. 4: The diffusion coefficient D versus x = Am/mo for iigiiter, heavier particies, and the totai. 
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FIG. 5: The velocity autocorrelation function for different Am. 
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FIG. 6: The velocity autocorrelation function ip{t) for lighter, heavier particles, and the total for Am = 0.7. 
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FIG. 8: a2(t) for lighter, heavier particles, and the total for Am = 0.7. 
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FIG. 10: a2 versus x = d^m/niQ. Stars represent results of computer simulation and the solid line is Eq. 
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FIG. 11: ln(Z)A/Ds) versus \n{mB/'mA)- Triangles are results of computer simulation and the solid line 
is the best fit line. 7 is the slope of this line. 



